ADAPTIVE, FAST AND OBLIVIOUS CONVOLUTION 
IN EVOLUTION EQUATIONS WITH MEMORY 
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Abstract. To approximate convolutions which occur in evolution equations with memory terms, 
a variable-stepsize algorithm is presented for which advancing steps requires only 0{N log N) op- 
erations and 0(logN) active memory, in place of 0(N'^) operations and 0{N) memory for a direct 
implementation. A basic feature of the fast algorithm is the reduction, via contour integral repre- 
sentations, to differential equations which are solved numerically with adaptive step sizes. Rather 
than the kernel itself, its Laplace transform is used in the algorithm. The algorithm is illustrated 
on three examples: a blow-up example originating from a Schrodinger equation with concentrated 
nonlinearity, chemical reactions with inhibited diffusion, and viscoelasticity with a fractional order 
constitutive law. 
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1. Introduction. We consider the problems of computing the convolution 
t 

f{t~T)g{T)dT, 0<t<T, (1.1) 

possibly with matrix- valued kernel / and vector- valued function g, and of solving 
evolution equations with memory containing such convolution integrals where g is 
not a function known in advance, but gir) depends on the solution at time r of the 
integral equation or integro-differential equation. In previous papers [16, 21] we have 
developed convolution algorithms that arc fast and oblivious: to approximate (1.1) on 
a grid t = nh (n = 0, 1, . . . , N) with constant step size h and Nh = T, the algorithm 
requires 

• O(A^logiV) operations, 

• 0(log A^) evaluations of the Laplace transform F ~ Cf, none of /, and 

• O(logA^) active memory. 

In the nth time step, g is evaluated at i„ = nh, but the history g{tj) for j < n is 
forgotten in this algorithm, and only logarithmically few linear combinations of the 
values of g are kept in memory. This is to be contrasted with the 0{N^) operations, 
0{N) evaluations of the kernel /, and 0(N) memory for a naive implementation of 
a quadrature formula for (1.1). Moreover, we note that in many applications the 
Laplace transform F (the transfer function), rather than the kernel / (the impulse 
response), is known a priori. A basic feature of the fast algorithm is the reduction, via 
contour integral representations, to differential equations of the form y' = Xy + g for 
suitable complex values of A, which are solved numerically. It is not necessary to solve 
these differential equations with constant time step h, as was done in [16, 21], but the 
step size may instead be adapted to the behavior of g. This observation opens the 
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way to an adaptive fast and oblivious convolution algorithm. Turning this simple idea 

into an cfBcicnt algorithm is, however, not as simple and the development of such an 
algorithm is precisely the topic of the present paper. The need to use adaptive time 
steps in solving evolutionary intcgro-differential equations in applications has been 
addressed at various places in the literature, e.g., by Adolfsson, Enclund & Larsson 
[4], Cao, Burrage & Abdullah [6], and Diethelm & Freed [8]. None of the adaptive 
algorithms proposed there, however, can make use of the convolution structure to 
reduce the 0{N'^) operation count and 0{N) memory requirements for N steps. The 
convolution algorithm proposed here works in the situation of a sectorial Laplace 
transform F: 

F is analytic in a sector | arg(s — cr)| < tt — (p with ip < ^n, 
and in this sector, \F{s)\ < M \s\~'' for some real M and u > 0. 

An equivalent condition is that / is analytic in a complex sector containing the posi- 
tive real half-axis t > 0, and is boimded by < C V''^^e'^*- in this sector. A typical 
example is the fractional-power kernel f{t) = t'^~^/T{i'), which has the Laplace trans- 
form F{s) = s"". An essential ingredient of the algorithm is the discretization of the 
inversion formula for the Laplace transform, given by 



with r a contour in the sector of analyticity oriented with increasing imaginary part 
and going to infinity with an acute angle to the negative real half-axis, so that e*'^ 
decays fast for growing |A| along F. We will choose the contour as a hyperbola. Since 
we cannot obtain a uniformly good approximation for all t G (0, T] with a single 
contour F, we use different hyperbolas F^ corresponding to geometrically growing 
intervals of uniform approximation, t G B^+^/i*] with an integer B > 1, e.g., 

i? = 5, and with a minimum step size h^,. The required number of contours is thus 
bounded by logg{T/hi:). This logarithm shows up in the complexity estimates in place 
of logg(Ar) for the fixed-stepsize algorithm. As in that case, it appears multiplied with 
the number of quadrature points on each hyperbola, which is 0(log -) for an accuracy 
s in the approximation of (1.3). In §2 we briefly review recent results from [14, 15] on 
the approximation of inverse Laplace transforms by discretized contour integrals. In §3 
we describe the fast and oblivious algorithm for computing convolutions with variable 
time steps. The algorithm is then illustrated on various problems where adaptive 
time steps are important: a blow-up problem for a nonlinear Abel integral equation 
resulting from a nonlinear Schrodinger equation with concentrated nonlinearity (§4), 
a fractional diffusion-reaction system from chemical reaction kinetics (§5), and visco- 
elasticity with a weakly singular memory kernel in the constitutive equations, under 
applied forces that are switched on and off (§6). 

2. Preparation: Numerical inversion of the Laplace transform. In the 

inversion formula (1.3) we choose F as the left branch of a hyperbola parameterized 



where /x > is a scale parameter, a is the shift in (1.2), and < a < tt/2 — ip. Thus, 
F is the left branch of the hyperbola with center at (/i, 0), foci at (0, 0), (2/i, 0), and 
with asymptotes forming angles ±(7r/2 -|- a) with the real axis, so that F remains 
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t>0 



(1.3) 



by 



M ^ F : X I— > 'y{x) = /i(l — sin(a -|- ix)) + a, 



(2.1) 



in the sector (1.2) of analyticity of F. After parameterizing (1.3), the function / 
is approximated by applying the truncated trapezoidal rule to the resulting integral 
along the real axis, i.e., 

1 r ^ 

fit) = 7r^ e'^ F{X) dA « y2 Wk e'^" F{\u), (2.2) 
with weights Wk and quadrature nodes Afc given by 

2m 

and r > a suitable step length parameter. Different choices of contours T and 
parameterizations have been studied for the numerical inversion of sectorial Laplace 
transforms in the last years. The choice of a hyperbola has been studied in [14, 15, 
17, 22, 9, 10, 28], and actually we follow here the approach in [14, 15]. The choice 
of r as a parabola has also been considered recently in [9, 10, 28]. Finally we refer 
to Talbot's method [26, 19, 27], which could also be used with the present algorithm; 
cf. [20, 16, 21, 13]. The good behavior of this quadrature formula to approximate 
(1.3) is due to the good properties of the trapezoidal rule when the integrand can 
be analytically extended to a horizontal strip around the real axis [24, 25]. We refer 
to [14, 15] for details and only give here the following error bound, which decays 
exponentially in the number of quadrature nodes. 

Theorem 2.1. [15] Suppose that the Laplace transform F satisfies the sectorial 
condition (1.2). For fixed T >0, A>1, 0<a< 7r/2 — and K > 1 there are positive 
numbers Ci, C2, C, c depending on a and A (C depends additionally on T unless a < 
in (1.2)^ such that the choice of parameters t = Ci/K and fj, = C2K/{Ato) yields a 
quadrature error in (2.2) bounded by 

\EKit)\<Ct''-' (e + e--^), 

uniformly for t in intervals [to,Ato] with arbitrary < to < T/A, where v is the 
exponent of (1.2) and e is the precision in the evaluations of the Laplace transform F 
and the elementary operations in (2.2). 

Hence, K — 0(log i) quadrature points are sufficient to obtain an accuracy 0(e) in 
the approximation of the contour integral. In practice, we choose a « ^(-1 — </?) and 
compute the values Ci and C2 following the optimization process described in [15]. 

3. The variable-step-size, fast and oblivious convolution algorithm. 

3.1. Local reduction to differential equations. We want to approximate 

u{t)^ f f{t-T)g{T)dT (3.1) 
Jo 

on a sequence of times < ti < ■ ■ ■ < tjq , spaced arbitrarily. For the moment we 
assume that 3 is a known function, though we will see later how to use the algorithm 
for solving integral and integro- differential equations. For a given we can insert 
the Laplace inversion formula in (3.1) and write 

^ " /(in - r) <?(r) dr = ^ " 2^ ^ e(*''-")^F(A) dA <?(t) dr . (3.2) 



The numerical inversion of the Laplace transform is performed very efficiently by the 
quadrature rule (2.2). In Section 2 wc have seen that the same contour T used in 
this quadrature can be used to approximate / at different values of t, ranging over 
intervals of the form [to,Ato], for a given ratio A > 1. Since in (3.2) we need to 
approximate /(f„ — t) for i„ — t S [0,t„], wc cannot use a imiquc contour T and we 
need to split the integral in (3.1) into several pieces. For suitable intermediate times 
< t~ < < tn, with {tn — t~)/{tn — t'^) < A, we select a suitable contour F for 
the time interval [tn — t~,tn — t'^] and approximate 

/* f{tn-T)g{T)dT = / ^/ e^*"-^^^F{X)dXg{T)dT 

/.t+ K 

« / Wke^*--^^^''F{Xk)g{T)dT 

k=-K 

= WkF{Xk)e^*"-'^^^'' / e(*+-")^'=5(r)dr 

k=-K 
K 

= J2 WkF{Xk)e^'''-'^'>^''y{t+,t-,Xk), (3.3) 

k=-K 

where y{t~^, t~,Xk) is the solution at t+ to the linear inhomogeneous ODE 

y' = Xky + g, yit-) = 0, -K<k<K. (3.4) 

We now approximate y{t'^ , Xk) by interpolating g linearly on each interval [tj, tj+i] 
for j = 0, . . . , — 1, and integrating exactly. (More elaborate integration methods 
could be used instead, cf. [7, 21], but for simplicity of presentation we will just consider 
this particular integration scheme.) We denote by g the interpolant of g and by 
y{t~^, t~,Xk) the resulting approximation to y{t~^, t~,Xk), i.e., y{t~^, t~,Xk) is the exact 
solution at t+ to 

y' = hy + 9, y{t-) = o, -K<k<K. (3.5) 

Thus, we approximate 

f{tn-T)g{T)dT^ Y WkF{Xk)e^'--'^^^'y{t+,t-,Xk). (3.6) 

k=-K 

3.2. Filling the mosaic. The key to the algorithm is the way the splitting 
times for the integral in (3.1) are selected for every t„ with 1 < n < N. This is 

done following the mosaic in the triangle {{t,T) : < r < i < T} shown in Fig. 3.1, 
where patches grow geometrically with increasing distance from the diagonal. For the 
moment, we fix a minimum size of the patches closest to the diagonal, corresponding 
to a minimum step size /i,. If along the vertical line at t„ joining with the diagonal 
value t„ we have L different patches of the tessellation, then we obtain the values tj 
and ior 1 < £ < L as the smallest and largest grid points, resp., within the £th 
patch along the vertical line at t„. All those £s are collected in an index set J. In 
case a patch does not contain any grid points its value i is not contained in J. The 
times tf = t depend on n, though for simplicity we omit this dependence in the 
notation. Each class of patches of the same size in the mosaic represents a distance 
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class to the diagonal in the mosaic, and thus corresponds to a different approximation 
interval and a different contour to perform the inversion of the Laplace transform and 
to a different set of 2K + 1 scalar differential equations. The approximation intervals 
for the values t„ — tf are of the form Ii = B^+^/i*], 2 < £ < L, so that 

the ratio A is B^. Since we consider a non-equidistant sequence of times tj, in this 
splitting there likely appear "gaps" in between the t^„j and tJ, which in Figure 3.1 
correspond to pairs of horizontal lines with any boundary of rn patches in between 
them. For example, at the time point ti5 = 3.45 we have L = 3, J = {1,3}, ^3 = 0, 
t+ = til = 2.11, = ti2 = 3.14 and if = hs = 3.24. 
We split (3.1) into 2| J| + 1 parts 

«(i„) = tiW(i„) + ^ u^'\u) + i^^'\tn) (3.7) 

where 

U^^\tn)= r fitn-T)g{T)dT (3.8) 

is computed using (3.3) and 

u^°\tn)= r f{tn-T)g{T)dT and u^^\tn)= f [tn - T)g{T) dr (3.9) 

Jt„-l Jit+m 
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correspond respectively to the step from tn-i to t„ near the diagonal and to the 
gaps between i^„j and tJ; sec Fig. 3.1. These parts are computed by "direct steps" 
explained in the next subsection. Thus, for ti5 = 3.45, ^(tis) is calculated using: 

• one "ode step", u^^\ from = to = to = tn. 

• one "direct step", vP\^ from = tn to = ti2- 

• one "ode step", u^^\ from = t\2 to if = ti^. 

• one "direct step", u^^\ from = tiz to tu. 

• one "direct step", u^"-*, from ii4 to ti^. 

Translating the above splitting from the picture into a formal procedure, we thus 
proceed as follows: given a minimum step size /i* and a basis B, for each tn we take 
L as small as possible so that we can represent [t„//i*] = 2 + X^fci^^-'^^"^ with 
hi e {1,2,..., B}, where \x\ denotes the smallest integer greater than x. is the 
largest and tJ is the smallest value in {tj : j = 0, . . . , n}, such that 

L L 

K X] bkB''-^ <tj <tj <kY^ bkB''-\ 
k=e+i k=e 

We remark that with this definition, there always exists some j G {0, . . . , n — 1} and 
some integer m > such that = tj and = tj+i- The tf are such that 

the integration limits of u^^\tn) fit into the approximation interval: tn — if G I^. 
In Figure 3.1 we show how the solutions to the linear ODEs y{t'[ ,tf^ , Xl.) fill the 
mosaic. In this example we have _B — 3 and the tj are given by the non-equidistant 
sequence of time points indicated both along the horizontal axis and the diagonal of 
the triangle. We can see here that, for instance, when times greater than t = 3.2 are 
reached, all the "past" from to 2.11 is stored in the solution values y(2.1, 0, A^^-*), 
represented in Figure 3.1 by the tall dark rectangle with basis [3.2, 4.1] for £ = 3, and 
in the adjacent incomplete white rectangle to the right for £ = 4. In our example 
the values y(2.1,0,A^"0 

are used 8 times to evaluate u{t) for t G [3.2,4.1]. The 
filling of the mosaic is done bottom up in the algorithm, advancing all the differential 
equations in every time step, so that the algorithm can forget all the past values of 
the function g, with the exception of those at tf, which are needed for the direct 

steps described below. In addition to the current solution values y{t, ,xf^) of the 
differential equations at time t = tn, also their values at a splitting point tf need to 
be stored until ff is increased at a later step. Actually the algorithms stores three 
copies of y{t'l ,tj , A^f^). Pseudocodes for the organization of the decomposition are 

given in the appendix. 

3.3. Direct steps. The gaps [i^„, , ^7] between the enclosed blocks in Figure 3.1 
are bridged using the values u^^\tn) whose computation we describe next. These 
direct steps compute 

^^'\tn) = I ' f{tn - T)g{T) dr = f - r)5(r) dr (3.10) 

for some j G {0, 1, . . . , N — 1}. On the interval [tj,tj+i] we approximate g{t) by a 
linear function: 

9{t) « m = 9j + ^^±i^(t - tj), hj+1 = tj+1 - tj, 

Hj+i 
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with Qj = g{tj), j = 0,1, . . . ,N. (Here again, the approach would extend to polyno- 
mials of higher degree.) Extending g{t) to [0, f„] we split (3.10) in two terms 



f{tn-T)-g{T)dT= r f{tn-T)g{T)dT- / " /(t„ " t) ^(t) dr 
Jtj Jtj Jtj+l 

= C-^[F ■ Cg{- + tj)]{tn - tj) - C-^[F ■ Cg{-+tj+i)\{tn - tj+i) 



h 



(in ij ) ^ 



^1 9j+l + ^2 7 



{tn - tj+l). 



where -Fi(s) = F{s)/s and F2{s) = F{s)/s'^. We approximate the inverse Laplace 
transforms 

Mt) = {jC-'F,){t), f2{t) = {C-'F2){t) 

at t = tn — tj and t = tn — tj+i using the numerical integration of the Laplace inversion 
formula along the integration contours corresponding to the approximation intervals 
/^i and 1^2 such that tn — tj+i £ I^^ and t„ — tj G li^ . The result of the direct step 
is then calculated forming linear combinations 



f{tn - t) g{T) dT = h{tn - tj) gj + hitn - tj) + ^ 



r .„ , 

9i+i - 9j 



(3.11) 



fl{tn - tj+l) gj+l - f2{tn - tj+l) 



hi 



This is also used for the terms closest to the diagonal, tj+i = tn, where we note in 

addition that /i(0) = /2(0) = 0. 

3.4. Complexity. Given an arbitrary sequence of time points < ti < . . . < 
tN = T with the minimum step size /i* = m.mj{tj+i — tj), the above algorithm 
computes 

t 

f{t-T)g{T)dT, for t = ti,...,tN, (3.12) 



(with g the piecewisc linear intcrpolant of g) up to an error e using 

• L = 0(log ^) hyperbolas with 

• 2K + 1 = 0(log i) quadrature points on each hyperbola. 

The algorithm thus requires {2K + 1)L evaluations of the Laplace transform F{s) 
at the quadrature points and solves {2K + 1)L differential equations y' = x'f^y + g- 
As the algorithm proceeds, only three solution values need to be stored for each of 
these differential equations. In addition, at most 2L values of g need to bo kept in 
memory for the direct steps. In total, the active memory requirements are 0{LK) = 
0(log ^ log i) vectors of the dimension of g. The total operation count is 0{NLK) = 
0{N log ^ log ^). For the variable-step-size algorithm we thus obtain the complexity 
characteristics as stated in the introduction for the fast and oblivious fixed-step-size 
algorithm, with log N now replaced by log ^ . 

3.5. Adaptivity based on controlling the interpolation error. There are 
two sources of error in the algorithm. The first one is the discretization of the contour 
integral, which is well controlled. The second one is the piecewise linear interpolation 
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of g by g. Ignoring the error from discretizing the contour integrals, the algorithm 
thus computes (3.12) instead of (1.1). We control the error in g, which is bounded by 

\\m - gm <^hl max \\g"{T)\\ for < t < 

*n-l5:''"Stn 

Given a tolerance Tol, we propose the new step-size hn+i according to the criterion 

Chl+.y^ = 0.8 -Tol, (3.13) 

where 7" = ||5'"(tn)|| with the quadratic interpolant 5 to 5 at tn-2,tn-i,tn- The 
constant C is chosen as C* ~ ^ \.f(t)\ dt. Additionally the step-size is restricted to 
fulfill l/2hn < hn+i < '2hn- The proposed step size hn+i is tested by 

'^/ln+l7n+l < Tol, (3.14) 

where the new value g{tn + hn+i) is used in the computation of 7^+1 • If this condition 
is satisfied, then hn+i is accepted and we set tn+i = tn + hn+i, else we repeat the 
test with a reduced step hn+i determined from (3.13) with j'^+i in place of 7^'. If 
necessary, this procedure is repeated until (3.14) is satisfied. 

In the following sections we give three examples to show the performance of 
the algorithm, using also other strategies for controlling the step size in the time 
integration of intcgro-differential equations. However, wc point out that the above 
fast algorithm is independent of the particular step size selection strategy. The step 
size control is just the way we generate the sequence of time points. The minimum step 
size need not bo specified a priori and at time t„ the future time points tn+2, ■ ■ ■ 
need not yet be known in the algorithm. 

For the three examples provided, we ran the algorithm with basis B = 5, which 
gives approximation intervals of the type [to, 25to]. 

4. A blow-up example originating from a Schrodinger equation with 
concentrated nonlinearity. Adami & Teta [1, 2] consider a nonlinear Schrodinger 
equation with nonlinearity concentrated at x = 0, 

j^ = -AV + 7|V'l'"V'-'5c.=o, x€R,t>0, (4.1) 

with initial data tjj{x,0) = ipo{x) for a; e M. The equation can be given a rigorous 

formulation as an integral equation. With Duhamel's principle the wave function 
tp{x, t) can be expressed as the sum of {U {t)il)o) {x), that is, the solution at position x 
and time t of the free Schrodinger equation with initial data ij^o, and of a convolution 
with Izp'^z, where z{t) := '0(0, t) for < r < /;. The function z is the solution of a 
nonlinear, complex Abel-type integral equation 

z{t) +l4 f rrj-^ Kr)\'''<r) dr = a{t) , t>0, (4.2) 

where a{t) = (U{t)ipo){0) is the solution at x = of the free Schrodinger equation. 
We present the results of numerical experiments in a situation where the solution is 
known to have finite-time blow-up. We choose a = 1, and 
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Fig. 4.1. Real and imaginary part (light and dark gray) and modulus (black) of the solution z 
for 7 = —2, —2.05, —2.06 — 2.5 and for different tolerances. 



which corresponds to a Gaussian wave-packet tpoix) = tt ^/^e ^^/^ as initial data. 

For four different values of 7, 7 = —2, —2.05, —2.06, —2.5 the evolution of the 
solution is shown in Fig 4.1 for different tolerances. The thick lines correspond to 
solutions obtained with a tolerance of 10~^. The other two lines correspond to tol- 
erances of 10~^ and 2 • 10"''. Whereas for 7 = —2, —2.5 one cannot distinguish the 
different tolerances in the plot, for 7 = —2.05, —2.06 it is clearly visible, that choosing 
a too low tolerance will produce a wrong result. 

In Fig. 4.2 the step size is plotted over t. Clearly the adaptivity pays off to resolve 
the blow-up. In case 7 = —2.5 we observe for a tolerance of 10"'' step-sizes ranging 
from 10"^ to 10"^. 

Fig. 4.3 displays the absolute error at the final time, obtained with tolerances 
1 • 10-3, 5 • lO"'', 2 • 10-'', 1 • lO-'', 5 • 10-^ 2 • 10-^ 1 • 10"^, 5 • 10"^ 2 • 10"^ 1 • 10"*^ 
versus the total number of steps. The error is measured against a reference solution 
obtained with tolerance 10"^. The numerical inversion of the Laplace transform is 
performed as explained in Section 2, with a = 0.8, d = 0.7, and K = 50, which gives 
Ci = 6.567 and C2 = 0.066 (sec Theorem 2.1). For larger tolerances, good results 
can bo obtained with a smaller K, say K = 25. Taking if = 40, we only get small 
oscillations in the stcpsize for the smallest tolerance, tol = 10~^, in Figure 4.2, and 
no visible changes for less stringent tolerances. 
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Fig. 4.3. Error versus the number of steps, at t = 10 for 7 : 
t = 3.15 for 7 = -2.06, and at t = 0.47 for 7 = -2.5 
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5. Chemical reaction kinetics with inhibited diffusion. We consider three 
molecular species A, B and C, reacting as 

A + B-^C 

C-^A + B 

C-^A + P, (5.1) 

P being the resulting product. The diffusion of each of the species is anomalous. So 
we obtain a reaction diffusion equation with a memory term. A model like this, with 
three species and three reactions was considered in [6]. However we have chosen to 
follow [29] and to associate a memory with the reaction term. Thus to model this 
process the following system of integro-differential equations is considered: 

ill =dtd^"[KAui - kiUiU2 + (^2 + fc3)u3) 

U2 =dtd^°' (KAu2 - kiUiU2 + ^2^3) (5-2) 

U3 =dtd^"{KAu3 + kiUiU2 - (^2 + kz)u3), 

where A = dxx is the ID Laplacian with periodic boundary condition on [—5, 5] and 
denotes the fractional integral of order < a < 1, given by the Riemann Liouville 
operator 

= ^)""'5(r) dr, (5.3) 

for < a < 1. Integrating in time system (5.2), we get the integro-differential 
equation 

ui(t) - ui(0) =dt°' [KAui{t) - kiui{t)u2{t) + (^2 + ks)ui{t)\ 

U2{t) - M2(0) =8^" [KAu2{t) - kiUi{t)U2{t) + k2U3{t)] (5.4) 
U3{t) - U3(0) =d^" [KAusit) + kiUi{t)U2{t) - (fc2 + k3)U3{t)] . 

In this situation we have in the convolution terms the weakly singular kernel 

f(t) = f"~^ /T{a), with Laplace transform F{s) = s^". Wc approximate the solution 
to (5.4) by using an adaptive strategy similar to the one explained in Section 3.5 but 
replacing criterion (3.13) by 

Chl+ij'^ = 0.8 -Tol, 

and the test (3.14) by 

Chl^,jUi<Tol, 

where 7^ = ||g'(i„)|l, with g the hnear interpolant of g at tn-i and tn- Our choice for 
the different parameters is K = 0.5, fci = 1, ^2 = 2 and ks = 3, and we integrate up 
to T = 30. We fix a = 0.5 and consider smoothed step- like functions as the initial 
data. 

Setting u = [ui,U2-Uo,]'^ , the 3x3 identity matrix, and following the notation 
introduced in Section 3.3 for the direct steps of the algorithm, the discrete equation 
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approximating (5.4) is 

f2{hn) 



= [fi{hn)-^^^){Kh^S + R)u"-^+ kih{hn)e^ur^u^-^ (5.5) 

1 5(u(r)) ^ 

r(a) Jo (* - T) 



where u denotes the piecewise linear interpolant of u at times to,ti, . . . ,tn and, for 
M nodes in the spatial discretization and v a column vector of length 3M, we define 

g(v) = {Kl3(g)S + R)v + fci (e V1V2) , e = [-1, -1, 1]'^, 

with S the second order finite difference approximation to d^x with periodic boundary 
conditions and R the 3M x 3M matrix 

ik2 + k3)lM 

R = \ fcz/M 

-{k2 + k3)lM 

where Im is the M x M identity matrix. 

Fig. 5.1 shows numerical results and performance characteristics of the algorithm. 
Here we inverted the Laplace transform taking a = 1, = 0.5, and K = 40 quadrature 

nodes on the hyperbolas, giving Ci = 6.036 and C2 = 0.0739. Again, less stringent 
accuracy requirements demand fewer quadrature nodes. 

6. Dynamic fractional order viscoelasticity. 

6.1. Model. The fractional order linear viscoelastic constitutive equation for 
the stress cr considered in [4, 3] reads 

(T = (To{t)-^ f f{t-T)ao{T)dT, (6.1) 

Jo 

with stress tensor ctq and strain tensor e given by 

ao{t) = 2ne{t) + Xtr{e{t))I ; e{t) = ^{Vu+ {Vuf) , (6.2) 

where ^ and A are the Lame constants and < 7 < 1 is a given parameter. We refer 

to [4] for more details about the model. 

The basic equations for the displacement field u are 

pu{x, t) — V ■ cro{u; x, t) 

+ 7/ f{t-T)V ■ cro{u\x,T) dr = for a; e f2; t > 
Jo 

u{x,Q) = uo{x) for a; e fi (6.3) 

u{x, 0) = vo{x) for a; g 
u{x,t) =0 ioY X ^Td] t>Q 
CTo{u;x,t) ■ n{x) = b{x,t) for a; G r^v; ^ > 0]. 
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Fig. 5.1. Left: Solutions for the three species A, B, C, at different times (lighter lines for 
larger times). Right: Step size versus time t, error versus tolerance, and error versus number of 
steps. All for T = 30, K = 0.5, fei = 1, fea = 2, and kz = 3. 



Td denotes the Dirichlet boundary of the O and Fiv the Neumann boundary, where 
the boundary force b is apphed. 

On the Sobolev space V := {v G H^{fl)'^ : v = on Yd} the variational formu- 
lation reads as follows: Find u{t) G V such that 

/ pu{x,t) ■ tp{x)dx + / €{u;x,t) : C€{ip{x))dx — / b{x,t) ■ 'tp{x)da{x) 
Jq Jn Jtn 

-1 j •^^^~'^)(/ ^{u]X,t) ■.Ce{il}{x))dx- j b(x,T) ■ip(x)da(x)^ dr (6.4) 

= , yxp ev 

u{0) = uo; u{0)=vo{x), 
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where the tensor product is given by 



2 

e{u) : Ce{tp) = ^ 2iJ,oeij{u)eij{tp) + Xoejj{u)eu{'ip) . 

Equation (6.4) is discretized in space nsing linear finite elements. The mesh is gener- 
ated using Triangle [23] and the assembly of the mass and stiffness matrices M and A 
and the boundary force vector b is done following [5] . In contrast to [5] we have chosen 
not to use Lagrange multipliers to enforcx; the Diriehlet data, but to incorporate the 
Dirichlet data directly. Thus (6.4) results in the abstract integro-differential equation 

Mu{t) + Au{t) - b{t) = 7 / f(t- t){Au{t) - 6(r)) 

Jo 

u{0) = uo ; u(0) = vq. 
The kernel / in (6.1) is given by 

f{t) = -j^E^{-n, 0<a<l, (6.5) 
where Ea denotes the Mittag-Leffler function of order a, defined by 

oo j 

The Laplace transform F oi f vs given by 

^ ^ l + s« 

6.2. Adaptive step size controL The discretization of the fractional order 
viscoelastic problem yields a Volterra integro-differential equation of second order of 
convolution type, 

Mu{t) + Au{t) = 7 / f{t- t){Au{t - 6(t)) dr + b{t) =: c{t). 
Jo 

This is equivalent to 

(■)^(-°. V)(:)-('0^ <") 

Applying the transformations u ^ it = M^^^u, = M~^/^w, A ^ A = 

M-V2^M-V2 and c ^ c = M'^^c, we get 

(^) = ( ' ) 
\ V J \ -Au + c J' 



In what follows we drop the "s. The time discretization is done using the Stormer- 
Verlet scheme, which is explicit and symmetric and has good properties for the part 
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u{t) = —Au{t) + b{t) (without the memory term). The Verlet scheme for the above 
equation reads 



V„+l/2 = Vn + ^{-AUn + Cn) 

Un+1 = Un + hVn+1/2 (6.7) 

h . . . 
Vn+1 = Vn+1/2 + -(-AUn+1 + C„+i), 

where c„ « c{tn) is computed using the adaptive convolution algorithm explained in 
Section 3. Note that w„ is already known before we evaluate the c„ and thus scheme 
is explicit. In order not to lose the good properties of the Stormer-Verlct scheme a 
special step-size control is used, following [11, 12] . For the integrating controller we 
fix an accuracy parameter e (which can roughly be viewed as the square root of a local 
error tolerance). Our step-size density function should control v and ii, therefore we 
take 

a{u,v,t)=a{u,v,t)-^^^={\\vf + \\ii\\l)~^^'' (6.8) 
= {{-Av + cf{-Av + c) + {-Au + cfA{-Au + c))"^^^ 

Assuming that A is symmetric, the partial derivatives of a are 

au{u, V, t) = 2{Au - cfAA 

aviu,v,t) = 2{Av - c)'^A 

at{u,v,t) = -2{Av - cfc-2{Au - cf^ Ac. 

With this choice the step-size becomes approximately proportional to 1/ -^HuH -I- c\\v\\ . 
We have to take the A norm of u so that and are in the same units. We 
use the integrating controller of [12], [11, (VIII. 3. 2)] and set 

1 / " \ 1 

G(u, V, t) = tV(t(u, V, tV -Au + c{t) \ = {2{Av - cf'c), (6.9) 

a{u,v,t) \ ^ j Aa 

where for an evaluation, c and c at tn are approximated by divided differences using 
c„,c„_i,c„_2 (set to zero for negative subscripts). 

Transforming back to "non-hat" quantities and again assuming that M and A 
are symmetric one obtains 

a{u, V, t) = {AM-'^v - cfM-'^{AM-'^v -c) + {Au - cf M-'^AM-'^{Au - c) (6.10) 

and 

G = -^{2{AM-'^v-cfM-'^c). (6.11) 

With an accuracy parameter e, and starting with c_i = c_2 = cq = bo, ^-1/2 = 
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l/a{uQ, Vo,to) — eG{uo, fo,to)/2, we compute for n = 0, . . . 

Zn+l/2 = + £G{Un, Vn, tn) 

hn+l/2 = 



^n+1/2 
tn+1 =^71 + hn+l/2 



Vn+1/2 =Vn + 



n+1/2 



{-AUn + Cn) 



Un+1 =Un + hn+l/2M ^Vn+1/2 

f2{hn+l/2^ 



Cn+1 = 7 



hn+l/2 



{AUn+1 - bn+l) + ( /l(/l„+l/2) - "+^/^-* j (^y^ _ 



h. 



n+l/2 



+ 



^ " - T)(^u(r) + b{T))dT^ + bn+l 



Vn+l = Vn+1/2 + 



T'n+1/2 



{-AUn+1 + C„+i), 



(6.12) 



where /i„+i/2 is the new step-size proposed by the integrating controller. 

6.3. Numerical example. In the example the domain Q. has the form of a 
cantilever as shown in Fig 6.1. The Dirichlet boundary Vjj - the left vertical boundary 

of n is indicated by squares. The time-dependent boundary force b is applied to the 
right vertical boundary FjVi of Q - indicated by circles in Fig. 6.1. On T^^ (the upper 
and lower part of the boundary of Q) homogenous Neumann boundary condition is 
assumed. 




0.5 1 1.5 



Fig. 6.1. Cantilever with finite element mesh. Tu is indicated by small squares and FjVj by 
circles. 



As the initial data u{x,0)^ v{x,Q) we take the lowest mode of the static semi- 
discrete problem corresponding to the first order equation (6.6), i.e.. the eigenvector 
of the matrix there corresponding to the smallest eigenvalue. The boundary force is 
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given by 



b(x t) = l 20eV((2*-5)«-i)(i,i)T for 2 < i < 3;a; e T^, 
\ else. 



In the numerical example the order of the Mittag-Leffler function is a = 1/2. We 
set the density p = I and 7 = 0.3. Youngs modulus and Poisson ratio are E = 200, 
v = 0.3. Equivalently the Lame constants are /x = 76.9, A = 115.4. For the numerical 
inversion of the Laplace transform we took here a = 0.8, d = 0.7 and K = 35 
quadrature nodes, giving Ci = 6.225 and C2 = 0.097. 

The evolution from t — to i = 6 of the horizontal component (x coordinate) of 
the displacement field it, the velocity it, and the acceleration it recorded at the upper 
left corner of the cantilever is shown in Fig. 6.2. At t = 2 when the boundary force is 
applied, an abrupt change in the velocity and strong oscillations in the acceleration 
is observed. Furthermore Fig. 6.2 shows the evolution of the step-size /'-„+i/2 for five 
different precision parameters e = 1 • 10"^4 • 10^^ 1 • 10"'', 2 • 10^4,5 • lO""*. The 
integrating controller reduces the step-size by roughly a factor of ten at t = 2. 

The relative error at t = 6 in energy norm + HwljAf is calculated with respect 
to a reference solution obtained with e = 1 ■ 10^^. The left plot in Fig 6.3 shows the 
relative error versus the number of steps. The different solutions were obtained for 

£ = 4 • lo-^ 5 • lo-^ 7 • 10-5, 1 . 10-^, 1.5 • 10-'', 2 • 10-*, 3 • lo-^ 4 • 10-^, 5 • 10-^. 
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number of steps number of steps 



Fig. 6.3. Error versus number of steps and error versus cpu time, for < = 6. 

The right plot in Fig 6.3 shows the cpu time in seconds versus the number of steps, 
illustrating the essentially linear growth. 
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Appendix A. Pseudocodes for the algorithm. We describe one step of the 
algorithm from to a given new time t = tn- 

For all = 1, . . . , L and k = —K, . . . ,K the ODEs (3.5) corresponding to the 

X], with initial time are advanced to the new time t or restarted, depending on 
whether the horizontal line at height t fits the current patch or enters a new patch 
(see Section 3.2 and Figure 3.1), i.e., we compute y^^\tj ,t, A^^^), or we set tj =t and 
y^^\tj ,t,\^^^) = 0, using the pseudocode 1. 

In the following pseudocodes Y(l) denotes a structure storing: 

• Y(l) .data^ y, the solution of the ODE corresponding to the £-patch, 

• Y(l) .tini <— tj, the initial time of the ODE in the ^-patch, 

• Y(l) .tcur •*— tn, the new time (new final time of the ODE), 

• Y(l).b ^ hi, the number of the current step in the corresponding ^-patch 
of the mosaic, along the vertical line from {tn,Q) to {tn,tn) (see Figure 3.1). 
This value ranges between 1 and B. In the example of Section 3.2 we have 
for tn = ^15, 63 = 3, 62 = 1 and bi = 3. 

• Y(l) .tmin <— t^^.^ = hminJ2k=e-\-i^kB''~^ , the baseline of the i?-patch of 
the mosaic, along the vertical line from {tn,0) to {tn,tn)- 

• Y(l) .tmax ^ tif}ax — hniin X]fe=£ ^kB^~^ , the top line of the ^-patch of the 
mosaic, along the vertical line from (t„,.0) to (f„,,i„,) (tmir'' = ^^mhi)- 



• Yil) .uh ^ hmini^ + ^k=o^^] upper bound of the approximation interval, 



Y(l) .Ib^ /i, 



"mm 




• Yd) .gini ^ s{tj) and Yd) .gcur<— g{tn), the values of the inhomogcncity 
at the initial time in the i'-patch and at the enrrent time. 
In the rare case that tn is exactly on one of the horizontal lines of the mosaic, the 
structure Yd) is copied to YA(1) before restarting for bookkeeping purposes. 



Algorithm 1 Advance and restart the scalar ODEs. odesol 
for / = 1 to L do 

if tn > Y{l).tmin + Y{l).b * B^'-^) * hmin then 
if tn > Y{l).tmax then 
if tn = Y.tmax then 

Y{1) = odesadvance(y(Z),i„,5„_i,c/„) ; 
Y{l).b^B ; 
YA{l)=Y{l) ; 
end if 

restart the ODE Y{1) c.f. Algorithm 2 ; 
else 

Y{1) = odesadvance(y(/),i„,g„_i,g„) ; 

while t > Y{l).tmin + Y{l).b * * hmin do 

Y{l).b = Y{l).b+l; 
end while 
end if 
else 

Y{1) = odesadvance(y(Z),i„,5„_i,5„) ; 
end if 
end for 



Algorithm 2 restart the ODE 

Y{l).data = Q*Y{l).data ; Y{l).b = 1 ; 

while tn > Y{l).tmAn + Y{l).b * * hmin do 

Y{l).b=Y{l).b+l ; 
end while 

Y [1) .tini ~ t; Y (1) .tcur = t; Y{l).gini = Qn', Y{l).gcur = Qn 'i 
while tn > Y(l).tmax do 

Y{l).tmin = Y{l).tmax ; 

Y{l).tm,ax = Y{l).tmin + * hmin ; 
end while 



Algorithm 3 Advancing the ODE (exponential Euler method) odesadvance 
for k = —K, . . . ,K do 

Y{i).datak = Y{(.).datak + {ex.p{dt * A^^^) - l)/x[^^ 

*{xi^^ * Y{e).datak + Qn-i + {9n - gn-i)/dt/x'j^^) + [gn-i - gnM^l 

end for 

Y{i).gcur = Qn ; Y{tj.tcur = t ; 



To fill the mosaic botton-up, Algorithm 4 is used. There copying the structure 
Yd) to YMd) and YAd) is done by checking if the distances to the diagonal t„ — 
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Y(l).tmin and tn — Y(.l).tmax At the approximation interval I(. 



Algorithm 4 update routine update 
for £ from L downto 1 do 

if tn-YM{e).tmin-hmin*YM{e)MB(-^-'^^ > Y{£).lb Sztn-YM.tini < Y{£).ub 
then 

YT{e) = YM{£) ; 
end if 

if tn > Y{£).tmax then 

YA{£) = Y{£) ; 

YM{1) ^Y{1); 
else if tn > Y{£).tmin + hmin * YM.b * B^^'^^ then 

YM{£) ^ Y{£) ; 
end if 

iftn-YM{£).tmm-hmimYM{£).h*B'^^-^) > Y {£) .lb k tn-Y M .tint < Y{£).ub 
then 

YT{£) = YM{£) ; 
end if 

if tn-YA{£).tmin-hmin*YA{e).b*B'^^-^^ > Y{£).lb ktn-YA.tini <Y{£).ub 
k YA{£).tini > YT{£).tini then 

YT{£) = YA(£) ; 
end if 
end for 



YT(1) is updated by checking if either of the distances to the diagonal correspond- 
ing to YM(1) or YA(1) fit the approximation intervals. 



Algorithm 5 update routine part 2 update 
idv = [ ] ; ito = \ \ 

if tn - YT{L).tcur > Y{L).lb k tn - YT{L).tini < Y{L).ub then 

idv = [L, idv] ; 
else 

YT{L).tini = -inf ; YT{L).tcur = -inf ; 
end if 

for £ from L —1 downto 1 do 

if tn - YT{£).tcur > Y{£).lb k tn - YT{£).tini < Y{£).ub k YT{£).tini > 
YT{£ + l).tini k YT{£).tcur > YT{£ + l).tcur then 

idv = [£, idv] ; 
else 

YT{£).tini = -inf ; YT{£).tcur = -inf ; 
end if 

if tn - tn-1 > Y{£).lb ktn- tn-1 < Y{£).ub then 

ito = £ ; 
end if 
end for 

idv = [ito idv] ; 



Algorithm 6 puts together the; ncccsary "direct" and "odes" steps. It uses the 
ode solutions YT and the vector idv from Algorithm 5. idv stores the orders £ of the 
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YT required at t = 



Algorithm 6 Fast convolution evaluation, c/. (3.7) 
out — ; 

if idv is not empty then 
if length(ic?w) > 2 then 

if tn-i > YT{idv{2)).tcur then 
out = out+ 

directstep(i, YT{idv{2)).tcur, YT{idv{2)).gcur, Qn-i) ; (3.11) 
end if 

for II = 2 : length{idv) — 1 do 

if YT{idv{ll)).tcur > YT{idv{ll)).tini then 
e = idv{ll) ; 

out = out + Ef=-if w^k^Fi>'k^) exp(t - YT{e).tcur) * YT{l).datak ; (3.6) 
end if 

if YT{idv{U)).tini > YT{idv{ll + l)).tcur then 

out = out + directstep(f, YT{idv{ll + l)).tcur, YT{idv{ll)).tini, 
YT{idv{ll + l)).gcur,YT{idv{ll)).gini) ; (3.11) 

end if 
end for 

if YT {idv {end)). tear > YT.{idv{end))tini then 

i = idv {end) ; 

out = out + Y.k=-K Wk^Fi^k^)^Mt - YT{l).tcur) * YT{e).datak ; (3.6) 
end if 
end if 
end if 
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